import os
import scipy.io
from datetime import datetime, timedelta
import pandas as pd
import numpy as np
from scipy.stats import norm
import shutil
import re
from scipy.special import comb

def isfield(py_dict, field_name):
    """
    Check if a field exists in a Python dictionary, similar to MATLAB's isfield function.
    
    Parameters:
    py_dict (dict): The dictionary to check.
    field_name (str): The field name to check for.
    
    Returns:
    bool: True if the field exists, False otherwise.
    """
    return field_name in py_dict
    
def addtodate(date_obj, amount, unit):
    """
    模拟 MATLAB 的 addtodate 函数，向给定的日期时间添加指定的年、月、日、小时、分钟或秒。
    
    参数:
    date_obj (datetime): 要修改的日期时间对象。
    amount (int): 要添加的量。
    unit (str): 时间单位，可以是 'years', 'months', 'days', 'hours', 'minutes', 'seconds'。
    
    返回:
    datetime: 修改后的日期时间对象。
    """
    if unit == 'years':
        # 添加年份
        new_year = date_obj.year + amount
        return date_obj.replace(year=new_year)
    elif unit == 'months':
        # 添加月份
        new_month = date_obj.month + amount
        new_year = date_obj.year + (new_month - 1) // 12
        new_month = (new_month - 1) % 12 + 1
        try:
            return date_obj.replace(year=new_year, month=new_month)
        except ValueError:
            # 如果日期无效（例如 2024-2-29），则返回上个月的最后一天
            return date_obj.replace(year=new_year, month=new_month, day=1) - timedelta(days=1)
    elif unit in ['days', 'hours', 'minutes', 'seconds']:
        # 添加天数、小时、分钟或秒
        delta = timedelta(**{unit: amount})
        return date_obj + delta
    else:
        raise ValueError("Unsupported unit. Use 'years', 'months', 'days', 'hours', 'minutes', or 'seconds'.")
def python_datetime_list(matlab_datetime):
    '''
    范围为python的浮点数数组;
    '''
    row, col = matlab_datetime.shape
    #print(matlab_datetime.shape)
    py_datetime = np.full([row,col], np.nan);
    for i in range(0, row):
        for j in range(0, col):
            #print(matlab_datetime[i,j])
            tmp = python_datetime(matlab_datetime[i,j])
            py_datetime[i, j] = tmp.timestamp()

    return py_datetime

def python_datetime(matlab_datetime):
    """
    matlab的时间系统转换为,python的datetime结构数据
    """
    # MATLAB 的起始日期是 0000 年 1 月 1 日
    #start_date = datetime(1, 1, 1)
    #dt_float   = start_date.timestamp();
    #delt = matlab_datetime - 366 - 1;
    #dt = datetime(1970, 1, 1, 0 , 0 , 0);
    #start_date4python = dt.timestamp();
    matlab_to_python_days = 719162 + 366 + 1
    # 将MATLAB的日期数转换为Python的日期数
    python_datenum4day = matlab_datetime - matlab_to_python_days
    # Python的起始日期是1970年1月1日
    start_date = datetime(1970, 1, 1)
    # 计算目标日期
    target_date = start_date + timedelta(days = python_datenum4day)
    return target_date
def matlab_datetime(date_python):
    """
    将python的时间类型转换为matlab的浮点数型;
    """
    # Python的起始日期是1970年1月1日，MATLAB的起始日期是0000年1月0日
    # 0000年1月0日到1970年1月1日的天数差
    days_since_0000 = 719162
    # 将Python的datetime对象转换为Unix时间戳（自1970年1月1日以来的秒数）
    unix_timestamp = (date_python - datetime(1970, 1, 1)).total_seconds()
    # 将秒数转换为天数，并加上MATLAB起始日期的偏移
    matlab_datenum = unix_timestamp / (24 * 60 * 60) + days_since_0000 + 366 + 1
    return matlab_datenum
def python_datenum(year, month, day, hour=0, miniute=0, second=0):
    """
    要求输入为np数组比较好;
    只能做一个变量的
    模拟 MATLAB 的 datenum 函数，将日期转换为自 0000 年 1 月 1 日以来的天数。
    参数:
    date_obj (datetime): 要转换的日期对象。
    返回:
    float: 自 0000 年 1 月 1 日以来的天数。
    """
    # MATLAB 的起始日期是 0000 年 1 月 1 日
    start_date = datetime(1, 1, 1)
    
    # 计算两个日期之间的天数差
    target_datetime = datetime(year, month, day, hour, miniute, second)
    delta           = target_datetime - start_date
    
    # MATLAB 的 datenum 返回的是浮点数，表示天数和时间的小数部分
    return delta.days + 366 + delta.seconds / (24 * 60 * 60) + 1;

def python_addtodate(date_fload, amount, unit):
    '''
    将matlab的代表时间的浮点数转换为matlab的代表的下个月的时间的浮点数
    '''
    date_obj = python_datetime(date_fload);
    new_date_obj = addtodate(date_obj, amount, unit);
    return matlab_datetime(new_date_obj);
def getmonthdatainfo(iy,im,id, type4date = '%Y-%m-%d'):
    SDATE     = datetime(iy, im, 1);
    EDATE     = addtodate(SDATE, 1, 'months');
    EDATE     = addtodate(EDATE,-1, 'days');
    STR_SDATE = datestr(SDATE,type4date);
    STR_EDATE = datestr(EDATE,type4date);
    return STR_SDATE, STR_EDATE
def getmonthdatainfo_matlab(iy,im,id):
    SDATE     = datetime(iy, im, 1);
    EDATE     = addtodate(SDATE, 1, 'months');
    EDATE     = addtodate(EDATE,-1, 'days');
    STR_SDATE = matlab_datetime(SDATE);
    STR_EDATE = matlab_datetime(EDATE);
    return STR_SDATE, STR_EDATE

def matlab_datestr(date_mat, format_str=None):
    date_obj = python_datetime(date_mat)
    STR_SDATE = datestr(date_obj, format_str)
    return STR_SDATE

def datestr(date_obj, format_str=None):
    """
    模拟 MATLAB 的 datestr 函数，将日期时间转换为字符串格式。
    
    参数:
    date_obj (datetime): 要转换的日期时间对象。
    format_str (str, optional): 日期时间格式字符串。如果未提供，使用默认格式。
    
    返回:
    str: 日期时间字符串。
    """
    if format_str is None:
        # 默认格式类似于 MATLAB 的默认 datestr 格式
        format_str = '%d-%b-%Y %H:%M:%S'
    
    # 将 datetime 对象转换为字符串
    return date_obj.strftime(format_str)

def getnextmonthfirstdate(iy,im):

    SDATE     = datetime(iy, im, 1);
    EDATE     = addtodate(SDATE, 1, 'months');
    jy        = np.nan;
    jm        = np.nan;
    jd        = np.nan;
    jH        = np.nan;
    jM        = np.nan;
    jS        = np.nan;
    [jy,jm,jd,jH,jM,jS] = datevec(EDATE);
    return jy, jm;

def update4mat(inFile, DDATA):
    SA = {'DDATA' : DDATA}
    if os.path.exists(inFile):
        ts = scipy.io.loadmat(inFile)
        for key, value in SA.items():
           ts[key] = value
        scipy.io.savemat(inFile, ts, do_compression = True)

def save4daydisp(inFile, DDATE, DDATA, DERRO, TDATE, TDATA):
    SA = {'DDATE' : DDATE,'DDATA' : DDATA,'DERRO' : DERRO, 'TDATE' : TDATE, 'TDATA' : TDATA}
    #print(inFile);
    #print(os.path.exists(inFile))
    if os.path.exists(inFile):
          #print('Exist')
          ts = scipy.io.loadmat(inFile)
          for key, value in SA.items():
              ts[key] = value
          scipy.io.savemat(inFile, ts, do_compression = True)
    else:
        #print('NotExist')
        scipy.io.savemat(inFile, SA, do_compression = True)
    return 1

def mkdir4file(iRoot):
    # 检查文件夹是否存在
    if not os.path.exists(iRoot):
    # 如果文件夹不存在，则创建它
      os.makedirs(iRoot)


def getfileinfo(infile):
    ifile = os.path.basename(infile)
    iroot = os.path.dirname(infile)
    return iroot,ifile

def splitfilename(isfile):
    parts = isfile.split('_');
    stnm = parts[0];
    eqid = parts[1];
    date = parts[2];
    return stnm,eqid,date

def datevec(date_input):
    if isinstance(date_input, str):
       date_input = datetime.strptime(date_input,'%Y-%m-%d %H:%M:%S')
    if isinstance(date_input, datetime):
       return [date_input.year, date_input.month, date_input.day, date_input.hour, date_input.minute, date_input.second];

def str2matlabdate_list(date_input, dbtype = '%Y-%m-%d %H:%M:%S'):
    N = len(date_input)
    date_ouput = np.full(N, np.nan)
    for i in range(N):
        date_ouput[i] = str2matlabdate(date_input[i], dbtype)
    return date_ouput
def str2matlabdate(date_input, dbtype = '%Y-%m-%d %H:%M:%S'):
    if isinstance(date_input, str):
       date_input = datetime.strptime(date_input, dbtype)
    if isinstance(date_input, datetime):
       return matlab_datetime(date_input)
def researchconditionfile4eq(instrument_id, iRoot):

    # 获取仪器ID和年月
    
    matching_files = []
    #year_month = '{:04}-{:02}-01-{:04}-{:02}'.format(iy,im,iy,im)
    #folder_path = '{}{:04}{:02}/'.format(iRoot, iy, im)
    folder_path = iRoot
    # 检查文件夹是否存在
    if not os.path.exists(folder_path):
        print(f"文件夹 {folder_path} 不存在。")
    else:
        # 遍历文件夹中的所有文件
        for filename in os.listdir(folder_path):
            # 检查文件名是否同时包含仪器ID和年月
            if instrument_id in filename:
                file_path = os.path.join(folder_path, filename)
                matching_files.append(file_path)
        
        if matching_files:
            print("找到以下符合条件的文件:")
            for file in matching_files:
                print(file)
        else:
            print("未找到符合条件的文件。")
    return matching_files
def search_files_by_keyword(folder_path, keyword):
    matched_files = [];
    for root, dirs, files in os.walk(folder_path):
        for file in files:
            if keyword.lower() in file.lower():
                full_path = os.path.join(root, file)
                matched_files.append(full_path)
    if matched_files:
       print("找到以下符合条件的文件:")
       for file in matched_files:
           print(file)
    else:
       print("未找到符合条件的文件。")
    return matched_files

def search_files_by_keyword_two(folder_path, keyword1, keyword2):
    matched_files = [];
    for root, dirs, files in os.walk(folder_path):
        for file in files:
            if keyword1.lower() in file.lower() and keyword2.lower() in file.lower():
                full_path = os.path.join(root, file)
                matched_files.append(full_path)
    if matched_files:
       print("找到以下符合条件的文件:")
       for file in matched_files:
           print(file)
    else:
       print("未找到符合条件的文件。")
    return matched_files


def append_vector(v1, v2):

    if v1.size == 0:
       return v2
    if v2.size == 0:
       return v1
    if v1.ndim == 1:
       v1 = v1.reshape(-1, 1) if v1.shape[0] > 1 else v1.reshape(1, -1)
    if v2.ndim == 1:
       v2 = v2.reshape(-1, 1) if v2.shape[0] > 1 else v2.reshape(1, -1)
    if v1.shape[0] != v2.shape[0] and v1.shape[1] != v2.shape[1]:
       raise ValueError('向量形状不兼容,无法合并')
    if v1.shape[0] == 1:
       result = np.hstack((v1, v2)) if v2.shape[0] == 1 else np.vstack((v1, v2))
    else:
       result = np.vstack((v1, v2)) if v2.shape[1] == 1 else np.hstack((v1, v2))
    return result



def getULbynormal(Y, alpha):
    # 处理缺失值，获取非空值
    ind = ~np.isnan(Y)
    M = Y[ind]
    
    # 估计正态分布的参数
    mu_hat, sigma_hat = norm.fit(M)
    
    # 计算置信区间的上下限
    lower_bound = norm.ppf(alpha / 2, mu_hat, sigma_hat)
    upper_bound = norm.ppf(1 - alpha / 2, mu_hat, sigma_hat)
    
    return lower_bound, upper_bound


def transferlongname2shortname(EQPOS, dictS=None):
    dictS = ['库车', '拜城', '沙雅', '精河', '呼图壁', '昌吉', '若羌', '和静', '尼勒克',
             '沙湾', '新源', '托克逊','善县','尉犁', '轮台','于田','克什米尔','洛甫','和田',
             '伽师', '皮山', '乌什','俄罗斯','塔城','且末','尼玛','临沧','墨江','宁蒗','漾濞',
             '盈江', '保山', '德宏', '缅甸', '印度','双柏',
             '那曲', '日喀则', '错那','墨脱','米林','丁青','门源',
             '玉树', '杂多', '玛多','海西','德令哈','茫崖',
             '九寨沟','夏河','积石山','马尔康','石渠',
             '阿拉善','泸定','泸县','巴东','遂宁','文县','剑河','宁强','威远',
             '青川','青白江','铜仁','武隆','秭归',
             '周口','应城','石柱','前郭','宁江','唐山','朝鲜',
             '灯塔','彰武','肃北','乐山','阿克陶','芦山']
    N     = len(dictS)
    M     = len(EQPOS)
    #print(M)
    arr   = np.array(EQPOS)
    for i in range(M):
        iEQPOS = EQPOS[i]
        for j in range(N):
            jdictS = dictS[j]
            result = np.char.find(iEQPOS, jdictS)
            if result != -1:
               EQPOS[i] = jdictS

    return EQPOS

def outputabnormaltimeserial(EQDATE, AB2D, AB2V, days4span):
    M = len(AB2D)
    
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    al_bz4d = [];
    al_bz4v = [];
    for i in range(M):
        iAB2D = AB2D[i]
        iAB2V = AB2V[i]
        for j in range(N):
            iEQDATE = EQDATE[j]
            if (iAB2D + days4span) > iEQDATE and iAB2D <= iEQDATE:
                bz = bz + 1
                al_bz4d.append(iAB2D)
                al_bz4d.append(iAB2D + days4span)
                al_bz4d.append(iAB2D + days4span + 1)
                
                al_bz4v.append(iAB2V)
                al_bz4v.append(iAB2V)
                al_bz4v.append(np.nan)
                break;
    #print(al_bz4d)
    #print(al_bz4v)
    return al_bz4d, al_bz4v

def outputcorrectedeq(EQDATE, AB2D, EQDIST, days4span):
    M = len(AB2D)
    
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    al_bz4d = [];
    al_bz4v = [];
    for i in range(N):
        iEQDATE = EQDATE[i]
        iEQDIST = EQDIST[i]
        for j in range(M):
            iAB2D = AB2D[j]
            if (iAB2D + days4span) >= iEQDATE and iEQDATE > iAB2D:
                al_bz4d.append(iEQDATE)
                al_bz4v.append(iEQDIST)
                break;
    #print(al_bz4d)
    #print(al_bz4v)
    return al_bz4d, al_bz4v

#绘制1条预测曲线;
#无论有没有地震都需要绘制;  
def outputpredictline(EQDATE, AB2D, AB2V, days4span):
    M = len(AB2D)
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    #只用判断最后一个异常点在预测时段是否有地震
    s_iAB2D = AB2D[M-1]
    a_iAB2V = AB2V[M-1]
    e_iAB2D = s_iAB2D + days4span
    al_bz4d = [];
    al_bz4v = [];
    al_bz4d.append(s_iAB2D)
    al_bz4d.append(e_iAB2D)

    al_bz4v.append(a_iAB2V)
    al_bz4v.append(a_iAB2V)
    return al_bz4d, al_bz4v

def copy4self(sfile, dfile):
    if os.path.exists(dfile):
        print(f"目标文件{dfile}已存在, 将覆盖.")
    shutil.copy(sfile,dfile)

def is_float(s):
    pattern = r'^[+-]?\d*\.?\d+(?:[eE][+-]?\d+)?$'
    return bool(re.match(pattern,s))

def is_int4one2zero(s):
    pattern = r'^[01]$'
    return bool(re.match(pattern,s))

def is_valid_date(date_str):
    pattern = r'^\d{4}-(0[1-9]|1[0-2])-(0[1-9]|[12]\d|3[01])$'
    #print(pattern)
    #print(date_str)
    #print(not re.match(pattern, date_str))
    if not re.match(pattern, date_str):
        return False
    try:
        datetime.strptime(date_str, '%Y-%m-%d')
        return True
    except ValueError:
        return False

#绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, abnormal_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_data
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;
def outputpredictdescripsion(DDATE,  DDATA,  DERRO, MDATE, MVALUE, TimeSpan, 
                             EQDATE, EQRANK, EQDIS, EQPOS,
                             yuzhi,  days4span,
                             sel4dis, sel4rank,
                             ouFile = None):
    #在进行预报效能评估的时候，地震目录要根据数据查询;
    abnormal_date,           abnormal_value             = outputabnormaldescripsion(MDATE, MVALUE, yuzhi)
    abnormal_corrected_date, abnormal_corrected_data    = outputpredict4abnormal2corrected(EQDATE, abnormal_date, abnormal_value, days4span)
    abnormal_wrong_date,     abnormal_wrong_data        = outputpredict4abnormal2wrong(EQDATE, abnormal_date, abnormal_value, days4span)
    earthquake_corrected_date, earthquake_corrected_rank, \
    earthquake_corrected_dis,  earthquake_corrected_pos = outputpredict4earthquake4corrected(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    earthquake_wrong_date, earthquake_wrong_rank, \
    earthquake_wrong_dis,  earthquake_wrong_pos         = outputpredict4earthquake4wrong(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    R                                                   = len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (MDATE[0,-1] - MDATE[0,0])
    successAlamCount                                    = len(earthquake_corrected_date)
    missAlamCount                                       = len(earthquake_wrong_date)
    R0                                                  = ywGetR0(successAlamCount, missAlamCount)
    #写文件;
    
    str_sdate = matlab_datestr(MDATE[0,0],'%Y%m%d')
    str_edate = matlab_datestr(MDATE[0,-1],'%Y%m%d')
    ybzysj    = days4span * len(abnormal_date)
    #print(ybzysj)
    ybzsj     = MDATE[0,-1] - MDATE[0,0]
    #print(ybzsj)
    line = "";
    line = "预测的数据时段:{}至{}\n".format(str_sdate, str_edate)
    line = "{}预测策略1, 震中距: {:.1f}km, 震级选择: {:.1f}, 地震数量: {:d}, 预报有效时长: {:.1f}d\n".format(line, sel4dis, sel4rank, len(EQDATE), days4span)
    line = "{}计算公式: R = 地震报对数({:d})/地震数({:d}) - 预报占用时间({:.0f}d)/预报总时间({:.0f}d)\n".format(line, len(earthquake_corrected_date), len(EQDATE), ybzysj, ybzsj)
    line = "{}预测策略2, 阈值:{:.4f}, R值:{:.3f}, R0值:{:.3f}\n".format(line, yuzhi, R, R0)
    if R > R0:
       line = "{}预报效能: 通过\n".format(line)
    else:
       line = "{}预报效能:未通过\n".format(line);
    line = "{}异常情况:  指标数: {:d}, 异常数: {:d}, 报准数: {:d}, 虚报数: {:d}\n".format(line, len(MDATE[0]), len(abnormal_date), len(abnormal_corrected_date), len(abnormal_wrong_date))
    line = "{}地震情况:  地震数: {:d}, 报对数: {:d}, 漏报数: {:d}\n".format(line, len(EQDATE), len(earthquake_corrected_date), len(earthquake_wrong_date))
    line = "{}报准地震列表:\n".format(line)
    for i in range(len(earthquake_corrected_date)):
        eq_date = matlab_datestr(earthquake_corrected_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_corrected_rank[i], earthquake_corrected_dis[i], earthquake_corrected_pos[i])
    line = "{}漏报地震列表:\n".format(line)
    for i in range(len(earthquake_wrong_date)):
        eq_date = matlab_datestr(earthquake_wrong_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_wrong_rank[i], earthquake_wrong_dis[i], earthquake_wrong_pos[i])
    abnormalsdate = matlab_datestr(abnormal_date[-1],'%Y%m%d')
    abnormaledate = matlab_datestr(abnormal_date[-1] + days4span,'%Y%m%d')
    line = "{}日期:{}的最新异常可能预示未来有相应地震发生\n".format(line, abnormalsdate)
    line = "{}预测截止时间为:{}".format(line, abnormaledate)
    print(line)
    if ouFile == None:
       ouFile ='/home/mw/temp/predict4des.txt'
    print(ouFile)
    with open(ouFile, 'w', encoding = "UTF-8") as file:
        file.writelines(line)
        
    return True 
#绘图后应该输出异常时间;
    #对预报的检验, 异常次数, 报准次数, 虚报次数
    #异常结构: abnormal_date, abnormal_value
    #报准结构: abnormal_corrected_date, abnormal_corrected_data
    #虚报结构: abnormal_wrong_date, abnormal_wrong_data
    #对预报的地震, 地震次数, 报准次数, 漏报次数
    #总的地震结构: earthquake_date, earthquake_rank, earthquake_pos, earthquake_dis
    #报准地震结构: earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_pos, earthquake_corrected_dis
    #报准地震结构: earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_pos, earthquake_wrong_dis
    #R值 =  len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (DDATE[-1] - DDATE[0])
    #R0值=  一个函数
    #需要写一段话;
def outputpredictdescripsion4speed(DDATE,  DDATA,  DERRO, MDATE, MVALUE, TimeSpan, 
                             EQDATE, EQRANK, EQDIS, EQPOS,
                             yuzhi,  days4span,
                             sel4dis, sel4rank,
                             ouFile = None):
    #在进行预报效能评估的时候，地震目录要根据数据查询;
    abnormal_date,           abnormal_value             = outputabnormaldescripsion(MDATE, MVALUE, yuzhi)
    abnormal_corrected_date, abnormal_corrected_data    = outputpredict4abnormal2corrected(EQDATE, abnormal_date, abnormal_value, days4span)
    abnormal_wrong_date,     abnormal_wrong_data        = outputpredict4abnormal2wrong(EQDATE, abnormal_date, abnormal_value, days4span)
    earthquake_corrected_date, earthquake_corrected_rank, \
    earthquake_corrected_dis,  earthquake_corrected_pos = outputpredict4earthquake4corrected(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    earthquake_wrong_date, earthquake_wrong_rank, \
    earthquake_wrong_dis,  earthquake_wrong_pos         = outputpredict4earthquake4wrong(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    R                                                   = len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (MDATE[0,-1] - MDATE[0,0])
    successAlamCount                                    = len(earthquake_corrected_date)
    missAlamCount                                       = len(earthquake_wrong_date)
    R0                                                  = ywGetR0(successAlamCount, missAlamCount)
    #写文件;
    
    str_sdate = matlab_datestr(MDATE[0,0],'%Y%m%d')
    str_edate = matlab_datestr(MDATE[0,-1],'%Y%m%d')
    ybzysj    = days4span * len(abnormal_date)
    #print(ybzysj)
    ybzsj     = MDATE[0,-1] - MDATE[0,0]
    #print(ybzsj)
    line = "";
    line = "预测的数据时段:{}至{}\n".format(str_sdate, str_edate)
    line = "{}预测策略1, 震中距: {:.1f}km, 震级选择: {:.1f}, 地震数量: {:d}, 预报有效时长: {:.1f}d\n".format(line, sel4dis, sel4rank, len(EQDATE), days4span)
    line = "{}计算公式: R = 地震报对数({:d})/地震数({:d}) - 预报占用时间({:.0f}d)/预报总时间({:.0f}d)\n".format(line, len(earthquake_corrected_date), len(EQDATE), ybzysj, ybzsj)
    line = "{}预测策略2, 阈值:{:.4f}, R值:{:.3f}, R0值:{:.3f}\n".format(line, yuzhi, R, R0)
    if R > R0:
       line = "{}预报效能: 通过\n".format(line)
    else:
       line = "{}预报效能:未通过\n".format(line);
    line = "{}异常情况:  指标数: {:d}, 异常数: {:d}, 报准数: {:d}, 虚报数: {:d}\n".format(line, len(MDATE[0]), len(abnormal_date), len(abnormal_corrected_date), len(abnormal_wrong_date))
    line = "{}地震情况:  地震数: {:d}, 报对数: {:d}, 漏报数: {:d}\n".format(line, len(EQDATE), len(earthquake_corrected_date), len(earthquake_wrong_date))
    line = "{}报准地震列表:\n".format(line)
    for i in range(len(earthquake_corrected_date)):
        eq_date = matlab_datestr(earthquake_corrected_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_corrected_rank[i], earthquake_corrected_dis[i], earthquake_corrected_pos[i])
    line = "{}漏报地震列表:\n".format(line)
    for i in range(len(earthquake_wrong_date)):
        eq_date = matlab_datestr(earthquake_wrong_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_wrong_rank[i], earthquake_wrong_dis[i], earthquake_wrong_pos[i])
    abnormalsdate = matlab_datestr(abnormal_date[-1],'%Y%m%d')
    abnormaledate = matlab_datestr(abnormal_date[-1] + days4span,'%Y%m%d')
    line = "{}日期:{}的最新异常可能预示未来有相应地震发生\n".format(line, abnormalsdate)
    line = "{}预测截止时间为:{}".format(line, abnormaledate)
    print(line)
    if ouFile == None:
       ouFile ='/home/mw/temp/predict4des.txt'
    print(ouFile)
    with open(ouFile, 'w', encoding = "UTF-8") as file:
        file.writelines(line)
        
    return True
def statisicabnormaltimespan(abnormal_date, days4span):
    ybzysj = 0
    if len(abnormal_date) > 0:
        N = len(abnormal_date)
        sdate  = abnormal_date[0]
        edate  = sdate + days4span
        isdate = sdate
        iedate = edate
        iybzysj= 0
        for i in range(N):
            cdate  = abnormal_date[i]
            isdate = cdate
            iedate = cdate + days4span
            #sdate edate 进行比较
            if isdate  >= sdate and iedate >= edate and isdate <= edate:
                iybzysj = iedate - sdate
                edate   = iedate
                print(iybzysj)
                if i == N - 1:
                   ybzysj  = ybzysj + iybzysj
            elif isdate > edate:
                sdate   = isdate
                edate   = iedate
                ybzysj  = ybzysj + iybzysj
                iybzysj = edate - sdate
                print(iybzysj)
             
    return  ybzysj


           
def outputpredictdescripsion4year4lo(DDATE,  DDATA,  MDATE, MVALUE, 
                             EQDATE, EQRANK, EQDIS, EQPOS,
                             yuzhi,  days4span,
                             sel4dis, sel4rank,
                             ouFile = None):
    #在进行预报效能评估的时候，地震目录要根据数据查询;
    abnormal_date,           abnormal_value             = outputabnormaldescripsion4lo(MDATE[0], MVALUE[0], yuzhi)
    print("异常时间")
    print(abnormal_date)
    print("异常时间")
    abnormal_corrected_date, abnormal_corrected_data    = outputpredict4abnormal2corrected(EQDATE, abnormal_date, abnormal_value, days4span)
    abnormal_wrong_date,     abnormal_wrong_data        = outputpredict4abnormal2wrong(EQDATE, abnormal_date, abnormal_value, days4span)
    earthquake_corrected_date, earthquake_corrected_rank, \
    earthquake_corrected_dis,  earthquake_corrected_pos = outputpredict4earthquake4corrected(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    earthquake_wrong_date, earthquake_wrong_rank, \
    earthquake_wrong_dis,  earthquake_wrong_pos         = outputpredict4earthquake4wrong(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    R                                                   = len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (MDATE[0,-1] - MDATE[0,0])
    successAlamCount                                    = len(earthquake_corrected_date)
    missAlamCount                                       = len(earthquake_wrong_date)
    R0                                                  = ywGetR0(successAlamCount, missAlamCount)
    #写文件;
    
    str_sdate = matlab_datestr(MDATE[0,0],'%Y%m%d')
    str_edate = matlab_datestr(MDATE[0,-1],'%Y%m%d')
    ybzysj    = days4span * len(abnormal_date) #这个是错误的;
    #print(ybzysj)
    ybzysj    = statisicabnormaltimespan(abnormal_date, days4span)
    ybzsj     = MDATE[0,-1] - MDATE[0,0]
    R         = len(earthquake_corrected_date) / len(EQDATE) - ybzysj / ybzsj
    #print(ybzsj)
    line = "";
    line = "预测的数据时段:{}至{}\n".format(str_sdate, str_edate)
    line = "{}预测策略1, 震中距: {:.1f}km, 震级选择: {:.1f}, 地震数量: {:d}, 预报有效时长: {:.1f}d\n".format(line, sel4dis, sel4rank, len(EQDATE), days4span)
    line = "{}计算公式: R = 地震报对数({:d})/地震数({:d}) - 预报占用时间({:.0f}d)/预报总时间({:.0f}d)\n".format(line, len(earthquake_corrected_date), len(EQDATE), ybzysj, ybzsj)
    line = "{}预测策略2, 阈值:{:.4f}, R值:{:.3f}, R0值:{:.3f}\n".format(line, yuzhi, R, R0)
    if R > R0:
       line = "{}预报效能: 通过\n".format(line)
    else:
       line = "{}预报效能:未通过\n".format(line);
    line = "{}异常情况:  指标数: {:d}, 异常数: {:d}, 报准数: {:d}, 虚报数: {:d}\n".format(line, len(MDATE[0]), len(abnormal_date), len(abnormal_corrected_date), len(abnormal_wrong_date))
    line = "{}地震情况:  地震数: {:d}, 报对数: {:d}, 漏报数: {:d}\n".format(line, len(EQDATE), len(earthquake_corrected_date), len(earthquake_wrong_date))
    line = "{}报准地震列表:\n".format(line)
    for i in range(len(earthquake_corrected_date)):
        eq_date = matlab_datestr(earthquake_corrected_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_corrected_rank[i], earthquake_corrected_dis[i], earthquake_corrected_pos[i])
    line = "{}漏报地震列表:\n".format(line)
    for i in range(len(earthquake_wrong_date)):
        eq_date = matlab_datestr(earthquake_wrong_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_wrong_rank[i], earthquake_wrong_dis[i], earthquake_wrong_pos[i])
    abnormalsdate = matlab_datestr(abnormal_date[-1],'%Y%m%d')
    abnormaledate = matlab_datestr(abnormal_date[-1] + days4span,'%Y%m%d')
    line = "{}日期:{}的最新异常可能预示未来有相应地震发生\n".format(line, abnormalsdate)
    line = "{}预测截止时间为:{}".format(line, abnormaledate)
    print(line)
    if ouFile == None:
       ouFile ='/home/mw/temp/predict4des.txt'
    print(ouFile)
    with open(ouFile, 'w', encoding = "UTF-8") as file:
        file.writelines(line)
        
    return True
def outputpredictdescripsion4year4up(DDATE,  DDATA,  MDATE, MVALUE, 
                             EQDATE, EQRANK, EQDIS, EQPOS,
                             yuzhi,  days4span,
                             sel4dis, sel4rank,
                             ouFile = None):
    #在进行预报效能评估的时候，地震目录要根据数据查询;
    abnormal_date,           abnormal_value             = outputabnormaldescripsion(MDATE[0], MVALUE[0], yuzhi)
    print("异常时间")
    print(abnormal_date)
    print("异常时间")
    abnormal_corrected_date, abnormal_corrected_data    = outputpredict4abnormal2corrected(EQDATE, abnormal_date, abnormal_value, days4span)
    abnormal_wrong_date,     abnormal_wrong_data        = outputpredict4abnormal2wrong(EQDATE, abnormal_date, abnormal_value, days4span)
    earthquake_corrected_date, earthquake_corrected_rank, \
    earthquake_corrected_dis,  earthquake_corrected_pos = outputpredict4earthquake4corrected(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    earthquake_wrong_date, earthquake_wrong_rank, \
    earthquake_wrong_dis,  earthquake_wrong_pos         = outputpredict4earthquake4wrong(EQDATE, EQRANK, EQDIS, EQPOS, abnormal_date, days4span)
    R                                                   = len(earthquake_corrected_date) / len(EQDATE) - days4span * len(abnormal_date) / (MDATE[0,-1] - MDATE[0,0])
    successAlamCount                                    = len(earthquake_corrected_date)
    missAlamCount                                       = len(earthquake_wrong_date)
    R0                                                  = ywGetR0(successAlamCount, missAlamCount)
    #写文件;
    
    str_sdate = matlab_datestr(MDATE[0,0],'%Y%m%d')
    str_edate = matlab_datestr(MDATE[0,-1],'%Y%m%d')
    ybzysj    = days4span * len(abnormal_date) #这个是错误的;
    #print(ybzysj)
    ybzysj    = statisicabnormaltimespan(abnormal_date, days4span)
    ybzsj     = MDATE[0,-1] - MDATE[0,0]
    R         = len(earthquake_corrected_date) / len(EQDATE) - ybzysj / ybzsj
    #print(ybzsj)
    line = "";
    line = "预测的数据时段:{}至{}\n".format(str_sdate, str_edate)
    line = "{}预测策略1, 震中距: {:.1f}km, 震级选择: {:.1f}, 地震数量: {:d}, 预报有效时长: {:.1f}d\n".format(line, sel4dis, sel4rank, len(EQDATE), days4span)
    line = "{}计算公式: R = 地震报对数({:d})/地震数({:d}) - 预报占用时间({:.0f}d)/预报总时间({:.0f}d)\n".format(line, len(earthquake_corrected_date), len(EQDATE), ybzysj, ybzsj)
    line = "{}预测策略2, 阈值:{:.4f}, R值:{:.3f}, R0值:{:.3f}\n".format(line, yuzhi, R, R0)
    if R > R0:
       line = "{}预报效能: 通过\n".format(line)
    else:
       line = "{}预报效能:未通过\n".format(line);
    line = "{}异常情况:  指标数: {:d}, 异常数: {:d}, 报准数: {:d}, 虚报数: {:d}\n".format(line, len(MDATE[0]), len(abnormal_date), len(abnormal_corrected_date), len(abnormal_wrong_date))
    line = "{}地震情况:  地震数: {:d}, 报对数: {:d}, 漏报数: {:d}\n".format(line, len(EQDATE), len(earthquake_corrected_date), len(earthquake_wrong_date))
    line = "{}报准地震列表:\n".format(line)
    for i in range(len(earthquake_corrected_date)):
        eq_date = matlab_datestr(earthquake_corrected_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_corrected_rank[i], earthquake_corrected_dis[i], earthquake_corrected_pos[i])
    line = "{}漏报地震列表:\n".format(line)
    for i in range(len(earthquake_wrong_date)):
        eq_date = matlab_datestr(earthquake_wrong_date[i], '%Y%m%d')
        line = "{}{:d} {} {:.1f} {:.2f}km {}\n".format(line, i + 1, eq_date, earthquake_wrong_rank[i], earthquake_wrong_dis[i], earthquake_wrong_pos[i])
    abnormalsdate = matlab_datestr(abnormal_date[-1],'%Y%m%d')
    abnormaledate = matlab_datestr(abnormal_date[-1] + days4span,'%Y%m%d')
    line = "{}日期:{}的最新异常可能预示未来有相应地震发生\n".format(line, abnormalsdate)
    line = "{}预测截止时间为:{}".format(line, abnormaledate)
    print(line)
    if ouFile == None:
       ouFile ='/home/mw/temp/predict4des.txt'
    print(ouFile)
    with open(ouFile, 'w', encoding = "UTF-8") as file:
        file.writelines(line)
        
    return True 
def ywGetR0(successAlamCount, missAlamCount):
    k = successAlamCount
    n = successAlamCount + missAlamCount
    alpha = 0.025
    minAlpha = float('inf')
    minAlpha_P = 0
    for P in [x * 0.001 for x in range(1001)]:
        sumAlpha = 0
        for i in range(k, n+1):
            sumAlpha += comb(n, i) * (P**i) * ((1 - P)**(n - 1))
        lessAlpha = abs(sumAlpha - alpha)
        if lessAlpha < minAlpha:
            minAlpha_P = P
            minAlpha   = lessAlpha
    R0 = k / n - minAlpha_P
    return R0
def outputabnormaldescripsion(MDATE, MVALUE, yuzhi):
    abnormal_date = []
    abnormal_value =[]
    indices = MVALUE >= yuzhi
    abnormal_value = MVALUE[indices]
    abnormal_date  = MDATE[indices]
    return abnormal_date, abnormal_value
def outputabnormaldescripsion4lo(MDATE, MVALUE, yuzhi):
    abnormal_date = []
    abnormal_value =[]
    indices = MVALUE < yuzhi
    abnormal_value = MVALUE[indices]
    abnormal_date  = MDATE[indices]
    return abnormal_date, abnormal_value
def outputpredict4abnormal2corrected(EQDATE, AB2D, AB2V, days4span):
    M = len(AB2D)
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    al_bz4d = [];
    al_bz4v = [];

    for i in range(M):
        iAB2D = AB2D[i]
        iAB2V = AB2V[i]
        for j in range(N):
            iEQDATE = EQDATE[j]
            if (iAB2D + days4span) > iEQDATE and iAB2D <= iEQDATE:
                bz = bz + 1
                al_bz4d.append(iAB2D)
                #al_bz4d.append(iAB2D + days4span)
                #al_bz4d.append(iAB2D + days4span + 1)
                al_bz4v.append(iAB2V)
                #al_bz4v.append(iAB2V)
                #al_bz4v.append(np.nan)
                break;
    abnormal_corrected_date = al_bz4d
    abnormal_corrected_data = al_bz4v
    return abnormal_corrected_date, abnormal_corrected_data
def outputpredict4abnormal2wrong(EQDATE, AB2D, AB2V, days4span):
    M = len(AB2D)
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    al_xb4d = [];
    al_xb4v = [];

    for i in range(M):
        iAB2D = AB2D[i]
        iAB2V = AB2V[i]
        isbz  = 0 # 一个异常只有虚报和报准两种情况 1 表示报准, 0表示虚报
        for j in range(N):
            iEQDATE = EQDATE[j]
            if (iAB2D + days4span) > iEQDATE and iAB2D <= iEQDATE:
                bz = bz + 1
                #al_xb4d.append(iAB2D)
                #al_bz4d.append(iAB2D + days4span)
                #al_bz4d.append(iAB2D + days4span + 1)
                #al_xb4v.append(iAB2V)
                #al_bz4v.append(iAB2V)
                #al_bz4v.append(np.nan)
                isbz = 1
                break;
        if isbz == 0:
           al_xb4d.append(iAB2D)
           al_xb4v.append(iAB2V)
        
    abnormal_wrong_date = al_xb4d
    abnormal_wrong_data = al_xb4v
    return abnormal_wrong_date, abnormal_wrong_data
def outputpredict4earthquake4corrected(EQDATE, EQRANK, EQDIS, EQPOS, AB2D, days4span):
    M = len(AB2D)
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    al_bz4d = [];
    al_bz4v = [];
    earthquake_corrected_date = []
    earthquake_corrected_rank = []
    earthquake_corrected_pos  = []
    earthquake_corrected_dis  = []
    for i in range(N):
        iEQDATE = EQDATE[i]
        iEQDIS  = EQDIS[i]
        iEQRANK = EQRANK[i]
        iEQPOS  = EQPOS[i]
        for j in range(M):
            iAB2D = AB2D[j]
            if (iAB2D + days4span) >= iEQDATE and iEQDATE > iAB2D:
                earthquake_corrected_date.append(iEQDATE)
                earthquake_corrected_dis.append(iEQDIS)
                earthquake_corrected_rank.append(iEQRANK)
                earthquake_corrected_pos.append(iEQPOS)
                break;
    #print(al_bz4d)
    #print(al_bz4v)
    return earthquake_corrected_date, earthquake_corrected_rank, earthquake_corrected_dis, earthquake_corrected_pos
def outputpredict4earthquake4wrong(EQDATE, EQRANK, EQDIS, EQPOS, AB2D, days4span):
    M = len(AB2D)
    #print(AB2D)
    N = len(EQDATE)
    #print(EQDATE)
    bz= 0
    earthquake_wrong_date = []
    earthquake_wrong_rank = []
    earthquake_wrong_pos  = []
    earthquake_wrong_dis  = []
    for i in range(N):
        iEQDATE = EQDATE[i]
        iEQDIS  = EQDIS[i]
        iEQRANK = EQRANK[i]
        iEQPOS  = EQPOS[i]
        isbz  = 0 # 一个异常只有虚报和报准两种情况 1 表示报准, 0表示漏报
        for j in range(M):
            iAB2D = AB2D[j]
            if (iAB2D + days4span) >= iEQDATE and iEQDATE > iAB2D:
                isbz = 1
                break;
        if isbz == 0:
           earthquake_wrong_date.append(iEQDATE)
           earthquake_wrong_dis.append(iEQDIS)
           earthquake_wrong_rank.append(iEQRANK)
           earthquake_wrong_pos.append(iEQPOS)

    return earthquake_wrong_date, earthquake_wrong_rank, earthquake_wrong_dis, earthquake_wrong_pos
#根据指标时段筛选地震目录
def selTimeSpanEqContent(MDATE, EQDATE, EQRANK, EQDIS, EQPOS):
    SDATE   = MDATE[0,0]
    EDATE   = MDATE[0,-1]
    #print(EQDATE)
    #print(SDATE)
    #print(EDATE)
    q    = np.where((EQDATE>=SDATE) & (EQDATE <= EDATE))
    print(q)
    oEQDATE = EQDATE[q]
    oEQRANK = EQRANK[q] 
    oEQDIS  = EQDIS[q]
    oEQPOS  = EQPOS[q]
    return oEQDATE, oEQRANK, oEQDIS, oEQPOS
    
# 示例用法
if __name__ == "__main__":
   date_input = '2024-01-01 00:00:00'
   #date_input = datetime(2024,1,1)
   date_obj = datevec(date_input)
   #print(date_obj)

















